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Quantifying the heart of darkness with GHALO - a 
multi-billion particle simulation of our galactic halo 
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ABSTRACT 

We perform a series of simulations of a Galactic mass dark matter halo at different 
resolutions, our largest uses over three billion particles and has a mass resolution of 
1000M o . We quantify the structural properties of the inner dark matter distribution 
and study how they depend on numerical resolution. We can measure the density 
profile to a distance of 120 pc (0.05% of i? v ir) where the logarithmic slope is -0.8 and 
-1.4 at (0.5% of i? v ir)- We propose a new two parameter fitting function that has a 
linearly varying logarithmic density gradient which fits the GHALO and VL2 density 
profiles extremely well. Convergence in the density profile and the halo shape scales as 
iV -1 / 3 , but the shape converges at a radius three times larger at which point the halo 
becomes more spherical due to numerical resolution. The six dimensional phase-space 
profile is dominated by the presence of the substructures and does not follow a power 
law, except in the smooth under-resolved inner few kpc. 

Key words: methods: N-body simulations - methods: numerical - dark matter - 
galaxies: haloes — galaxies: clusters: general 



1 INTRODUCTION 

Over twenty five years ago the theoretical framework for 
the evolution of a cold dark matter (CDM) dominated uni- 
verse was established i|Peebleslll983 ). The hierarchical and 
violent growt h of structure in this model begins at a scale 
of 10 _6 M Q (|Diemand et al.ll2005l ) until the most massive 
clusters of galaxies form that are many orders of magnitude 
more massive. The assumption that the dark matter is cold 
remains to be verified, yet numerical simulations that fol- 
low the hierarchical formation of CDM haloes have given 
several fundamental and robust predictions for the struc- 
tural and substructure properti es of the dark matter distri- 
bution within virial i sed haloes dDubinski fc Carlberdll99ll: 
Navarro et al.lll99rj; lAvila-Reese et a.1.1 Il999l : iBullock etail 



20011 : iFukushige fc Makinoll200ll ). These results are widely 



used to compare with observational data and to assist com- 
parisons with analytic models. 

The first CDM halo simulated with enough resolution t o 
resolve substructure used 10 6 particles (jMoore et al.l ll998). 
resolvi ng the density profile to about one percent of the virial 
radius (IGhigna et al.lll998l; iMoore et al.l l 19991 : iPower et all 
l2003l ; lNavarro et alj200l : lDiemand et al.ll2004T ). Whilst such 
simulations find numerous substructures in the outer halo, 



they find few or none within the inner 20% of _R v ir and no 
obvious structure i n phase-space in the central halo regions 
jMoore et alj|200lt) . Advances in algorithms and supercom- 
puting power have recently allowed us to increase this reso- 
lution by over two or ders of magnitude wi th the Via Lactea 
II (VL2) simulation (|Diemand et alj|2008j ). 

There are several reasons why we wish to carry out fur- 
ther studies at a higher resolution: (i) There are many old 
and forthcoming observational tests that constrain the struc- 
ture of dark matter haloes on scales well within 0.0017? v ir. 
These include high resolution rotation curve data and the 
kinematics of stars at the centres of dwarf galaxies. Future 
proper motions of these inner stars with GAIA or SIM will 
provide even tighter constraints. The close binary nuclei in 
galaxies such as VCC128 constra ins the dynamics on even 
smaller scales l|Goerdt et al.|[20"08h . (ii) As large surveys have 
pushed the surface brightness limits and detection efficien- 
cies, many extremely faint satellite galaxies have been found 
orbiting the Milky Way. The completeness of current sur- 
veys is debated, and it has been argued that many hundreds 
of additional system s may be found in the coming years 
ijTolIerud et alj|2008h . Simulations that can resolve and fol- 
low the survival of substructure within 10% of _R v j r are nec- 
essary to compare with these data, (iii) Dark matter de- 
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Figure 1. The density of dark matter within the inner 200 kpc of 
GHALO2. There are about 100,000 subhaloes that orbit within 
the virial radius. Each bright spot in this image is an individ- 
ual, bound, dark matter subhalo made up of many thousands of 
particles (there are far more particles than pixels here). 

tection, either directly on Earth or indirectly via detection 
of annihilation relics, is the ultimate way to determine its 
nature. These experiments rely on accurate predictions for 
the phase-space structure of dark matter at the position of 
the Earth's orbit and the abundance and inner properties 
of substructure throughout the Galactic halo, (iv) Under- 
standing the equilibrium structure resulting from violent re- 
laxation is the ultimate challenge for galactic dynamicists. 
There is no compelling theory that can explain universal 
densi ty and phase-space density profiles (|Tavlor fc Navarro] 
2001), or correlations such as betw een the local density pr o- 
file and the anisotropy parameter (jrlansen fc Sta del 2006). 

Given this motivation, we have carried out a sequence 
of simulations of a single Galactic mass dark matter halo, 
which at our highest resolution contains over a billion par- 
ticles within its virial radius. In this letter we report on its 
inner structure and convergence properties. 



2 THE SIMULATIONS 

Our initial condition s are based upon the WMAP3+SDSS 
l|Spergel et all 120071 : lYao et all 120061 ) cosmological model 
with erg = 0.742, Q M = 0.237, fi A = 0.763, h = 0.735, n = 
0.951. The galaxy sized, 10 12 M Q , i? v ir = 240 kpc, halo was 
selected from a cosmological cube of 40 Mpc on a side. 
This simulation had 488 3 particles (simulation GHALO5) 
in which three further nested spatial refinements by a factor 
of 3 (GHAL04,3,2) were placed such that the Lagrangian 
region of about 3i? v ir of the halo at z = was covered by 
2.1 x 10 9 high resolution particles in the initial condition. 
The final effective resolution of GHALO2 is 13176 3 result- 
ing in a particle mass of 1000 Mq and a total of 3.1 x 10 9 



particles and 1.3 x 10 9 particles within -R200 = 347 kpc. This 
allows us to capture all substructures out to more than 2i?200 
at the highest resolution. A further refinement GHALOi (in 
progress) will resolve the phase-space structure at the posi- 
tion of the sun more sharply for future recoil dark matter 
detection experiments. 

Creating these initial conditions was a significant chal- 
lenge an d we had to pa r alleliz e the GRAFIC1 and GRAFIC2 
codes of iBertschingeil (|200ll ) whereby the GRAFIC2 code 
was completely rewritten in C and MPI, and checked 
for near machine precision agreement with the original 
GRAFIC2. The new parallel GRAFIC1&2 codes can be ob- 
tained from the authors. Generation of the initial condition 
took 10 hours on 500 CPUs. We found that the original 
GRAFIC2 code had a bug in which the power spectrum 
used for the refinements was effectively that of the bary- 
onic component. Although this has affected many previous 
simulations (not GHALO, nor VL2), tests show that the 
conclusions of these studies are not compromised. 

The GHALO2 simulation was run at the Barcelona Su- 
percomputer Center on 1000 CPUs of Marenostrum using 
a total of 2 million CPU hours. Several significant improve- 
ments to the gravity code PKDGRAV2 made this calcula- 
tion possible including much better parallel computing effi- 
ciency and SIMD vector processing. P KDGRAV2 use s a fast 
multipole method (FMM) similar to iDehnenl (|2000l , |2002| ) 
but using a 5th-order reduced expansion for faster and more 
accurate force calculation in parallel, and a multipole based 
Ewald summation technique for periodic boundary condi- 
tions (|Stadelll200il ). It uses adaptive individual time-steps 
for parti cles based on a ne w estimator of the local dynami- 
cal time l|Zemp et alj|2007l ). The opening angle in the grav- 
ity tree and the accuracy parameter in the dynamical time- 
stepping is O = 0.55 and rj — 0.03 before z = 2, and then 
increased to 0.7 and 0.06 respectively. We make several com- 
parisons to the VL2 simulation which was also run with the 
FMM version of PKDGRAV2, but whose initial conditions 
were selected and generated independently using somewhat 
different methods. The VL2 halo has a mass of 2 x 1O 12 M 
and used a particle mass of 4000 Mq. The spline softening 
lengths for GHALO2, VL2, GHAL0 3 ,4, 5 are 61, 40, 182, 546, 
and 1639 pc, respectively (for GHALO these are set to 1/50 
of the mean inter-particle separation). 



3 THE INNER HALO STRUCTURE 

3.1 The dark matter density profile 

We apply a logarithmic binning to determine the radial 
density profile for the various simulations which are shown 
in Figure 2. The convergence radius of the density profile 
for the lower resolution realizations (GHAL03,4,5) can be 
clearly seen and are shown by the tick marks. These scale 
roughly as expected with r conv <x A?"" 1 / 3 , and we extrapolate 
this to conclude that the convergence radius of GHALO2 is 
around 120 pc. The inner slope of GHALO2 is -0.8 at 120 pc 
= 0.05% of Tivir and -1.4 at 2 kpc where the first subhalos 
become visible. Also shown is the power-law slope as a func- 
tion of log(r), which exhibits a similar linear functional form 
for both haloes with no rescaling. Based on this observation 
we propose a new functional form for the fitting function of 
the density profile, 
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A 2 (xlCT 4 ) p [lO 6 M fcpc -3 ] R[kpc] 3rd parameter 



Fitting Function 


Hernquist (a, /3, 7) (p s ,R s ) 


GH 2 


VL2 


GH 2 


VL2 


GH 2 


VL2 


GH 2 


VL2 


NFW 


(1,3,1) 


3.5 


6.6 


2.32 


4.24 


14.1 


13.9 






Dehnen-McLaughlin 


( 4 /9, 31 /9,79) 


1.6 


0.70 


0.273 


0.591 


42.6 


36.7 






S&M-profile {p ,Rx) 




0.93 


0.41 


5050 


11000 


2.20 


1.88 






Generalized NFW 


(1,3,7) 


3.0 


2.7 


1.78 


1.87 


16.2 


20.9 


1.04 


1.13 


Dehnen-McLaughlin 


(( 4 - 2 /3o)/ 9) ( 31 - 2 A))/ 9| (7+io/3 y g ) 


1.3 


0.68 


0.466 


0.522 


32.0 


39.1 


-0.0531 


0.0129 


Prugniel-Simien (p',R e ,a) 




1.5 


0.94 


14.0 


19.5 


59.6 


92.4 


0.376 


0.328 


Einasto (p_ 2 , ii!_ 2 , ») 




1.0 


0.45 


0.685 


0.991 


26.8 


28.9 


0.155 


0.142 


S&M-profile (p ,R x ,X) 




0.92 


0.41 


4710 


11200 


2.47 


1.82 


0.102 


0.100 



Table 1. Fitting parameters and A 2 for each of the 2 and 3-parameter models for both GHAL0 2 and VL2 simulations. Here A 2 = 
2^ m (ln(pi) — ln(pMODEL(''i))) 2 /("^ — 3) where p; are the density values in logarithmically spaced radial bins at . We fit from the resolved 
radiu s to 15% of -R v j r at which point substructure begins to cause significant fluctuations in the profile. Consistent with (Diemand ct al. 
IgOOSj) we obtain a generalized NFW with (p s ,R a ,^/) = (1.05,28.0,1.23) (units as above) for VL2 by fitting from 360 pc to R vil , with the 
best fit profile being Prugniel-Simien over this range, (p',R e ,ct) = (18.3, 113,0.308). 



p(r) = poe v v ' x " 



(1) 



which we term the S&M-profile (Stadel & Moore in prepa- 
ration). It is linear in this plot down to a scale R\ be- 
yond which it approaches the central maximum density 
po as r — > 0. We also note that if one makes a plot of 
din p/dln(l + r/R\) vs. ln(l-r-r/7?,\), then this profile forms 
an exact straight line with slope — 2A. 

Table 1 lists the best fitting parameters for several func- 



profiles ( 


Herat 


uist 199d; Zhaol 19961). the Einasto profile 


( Einasto 


1969; 


Navarro et al.ll2004l) 



p(r) = p_ 2 exp(-7 a [(rAR_ 2 r-l]), 
and the iPrugniel fc Simienl (|l997l ) profile 
p{r) = p'(r/R e )- p <* exp(-Mr/i^ 2 ) Q ), 



(2) 



(3) 



where p a = 1 - 0.6097a + 0.05463a 2 an d b a = 2 / a - % + 



0.009876a (for a < 2, see iMerritt et ail (|2006|)) such that 
when projected one obtains a Sersic profile (Sersic 1963, 
1 19681 ). 

The residuals shown in Figure 2 show that the S&M- 
profile provides a slightly better fit than all the models for 
the inner, more consistent, part of the profile. Furthermore, 
it is the only 3-parameter model where the 3rd parameter 
has a consistent value for the two different simulations. For 
this reason we also list this model as a possible 2-parameter 
model, fixing A = 0.1. The Einasto profile also provides an 
excellent fit to the density profiles of the two simulations. 



3.2 Convergence of Halo Shape 

The convergence o f the shape parameters (see also 
lAllgood et all (|2006l )) for GHALO in Figure 3 show that 
it is highly prolate over all resolved regions with b/a = c/a 
~ 0.5. At the halo centres the shape diverges quickly to a 
more spherical configuration. This is likely due to the orbital 
distribution being modified by the effects of resolution and 
softening. In this region the velocity distribution function is 
also strongly affected. 

We estimate the convergence in the shape to be achieved 
at 0.3, 0.6, 2, 15 kpc for GHAL0 2j 3,4,5 respectively, a radius 
that is about 3 times the inferred convergence radius of the 
density profile but also scaling as N' 1 ^ 3 . The fact that the 
variation in shape has little impact on the density profile 



can be understood by comparing the density profile taken 
in a 15 degree cone about the major, a, axis and the minor, 
c, axis (|jing fc Sutdl2002l ). The A 2 for the fits to the various 
density profiles remains roughly consistent between the two 
axial density profiles, although the best fit parameters vary. 
Due to the prolate shape the density profile parameters for 
the short axis are similar to the ones presented in Table 1. 



3.3 Phase-Space Density Profile 



It ha s been pointed out dTavlor fc Navarrd l200ll : iDehnenl 
120051 ; iDehnen fc McLaughlin! 120051 ) that the phase-space 



density (PSD) proxy, pa~ vs R is a power-law for CDM 
haloes, and several new fitting functions for the density pro- 
file have been proposed using this fact as a starting point 
such as the Dehnen-McLaughlin models. When averaged in 
shells, p(2na 2 )~ 3 ^ 2 is remarkably well fit by a power-law 
with slope of —1.84 as shown in Figure 4. However, it is in- 
teresting to compare this spherically averaged estimate with 
the true 6-dimensional PSD. 

The code EnBiD (ISharma fc Stcinmctz 200(3) has im- 
proved on earlier work bv lAscasibar fc Binnevl (|2005l ) in cal- 
culating better estimates of the 6-dimensional phase-space 
volume occupied by each particle and hence the PSD. Tak- 
ing the mean EnBiD PSD in logarithmic shells we see that 
the closest subhalo at 1.8 kpc stands out prominently and 
subhalos at larger radii begin to dominate the mean. Us- 
ing a method based o n a 6-dimensional Voronoi tessella- 
tion lArad et all (|2004l ) also showed that the subhalos form 
a dominant contribution to the phase-space density. This 
feature of using the EnBiD PSD can be turned to great 
advantage in identifying subhalos and other substructures 
such as phase-space streams. However, removing the effect 
of subhalos with /EnBiD > 100 M0kpc -3 (km/s)~ 3 from the 
mean, we extend the mean background PSD out to much 
larger radii as shown in Figure 4. By removing streams, 
with /EnBiD > O.4M0kpc -3 (km/s) ~ 3 , we can extend this 
to at least 40 kpc. 

We find that the true radial PSD profile estimated with 
EnBiD does not follow such a perfect power law and shows 
a steeper slope (roughly —2) than the per -3 estimator. The 



EnBiD mean estimate and p{2-kg 



-3/2 



are in agreement 



from about 0.2 to 2 kpc but the meaning of the power-law 
behaviour of p(27rer 2 ) -3 / 2 is unclear given that inside of 0.2 
kpc it is under-resolved and outside of 2 kpc a large con- 
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Figure 2. The upper panel shows the density profile of GHALO2 
and its lower resolution realizations as well as the density pro- 
file of the VL2 simulation in magenta. The convergence radius 
at each step in resolution is easily seen (indicated by the tick 
marks). The lower panel shows the residuals of the GHALO2 
simulation with respect to 2-parameter fitting functions: NFW 
(blue) and Dehnen-McLaughlin (green); as well as 3- parameter 
fitting functions: S&M-profile (black), Einasto (red), Generalized 
NFW (cyan), Dehnen-McLaughlin (magenta), Prugniel-Simien 
(yellow) . 
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Figure 3. Shape parameters for GHAL02,3,4,5- 
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Figure 4. The phase-space density profile of the main halo, mea- 
sured in several different ways is shown. The solid green line 
shows the traditional p(2ira 2 )~ 3 / 2 averaged in shells. The solid 
black and dashed black curves show s the mean and median En - 
BiD phase-space density estimator l|Sharma fc Steinmetj I2006F) 
for the particles in logarithmic shells extending out to 40 kpc. 
The blue and cyan curves show the mean EnBiD phase-space 
density profile, but where the contribution from subhalos (blue) 
and subhalos+streams (cyan) has been excluded. Despite the ef- 
fects of substructure the per -3 profile is remarkably well fit by a 
power-law with slope —1.84 shown as the dotted black line. 



tribution comes from the substructure. A further concern is 
the considerable variation of po~ 3 about a spherical shell 
of the prolate inner halo, which makes it remarkable that 
we obtain the same pow er-law slope as originally found by 
jTavlor fc Navarroll200ll ) despite the averaging that is tak- 
ing place. This also explains the good performance of the 
Dehnen-McLaughlin 2 and 3-parameter models at fitting the 
density profile. 

From about 2 to 40 kpc the po~ i estimator is somewhat 
enhanced due to the presence of substructure, while inside of 
0.1 kpc the EnBiD-mean continues to resolve the power-law 
behaviour of the profile. 



4 CONCLUSIONS 

The GHALO2 simulation has achieved an unprecedented 
spatial and mass resolution within a CDM halo, resolving 
thousands of subhalos within a radius corresponding to the 
galactic disk and a rich phase-space structure of streams be- 
yond a radius of ~ 8 kpc. Whilst there are more detailed 
analyses of this simulation in progress, we have reported 
here on the global inner properties of density and phase- 
space density profiles and halo shape. Using a sequence of 
simulations of the same halo at difference resolutions, from 
10 5 - 10 9 particles, we confirm that the convergence radii for 
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the density profile and shape scales as N vil ' . The logarith- 
mic slope of the radial density profile is close to a power law, 
gradually turning over to a slope of —0.8 at our innermost 
resolved region (0.05% of -R v i r ). We have proposed a new two 
parameter fitting function that has a linearly varying loga- 
rithmic gradient which provides the best fit to the inner part 
of the GHALO and VL2 ha loes. A larger sample of haloes, 
such as lBullock et all (|200lT > and lMaccio et all l|2007l ). would 
be required to determine if this functional form provides a 
universal fit. We find that the convergence radius of the 
density is a factor of three smaller than the convergence of 
halo shape. GHALO is prolate, yet becomes spherical within 
a region where orbits are most likely innacurately followed 
due to the effects of finite particle number, relaxation and 
softening. 

All functional forms fit to density profiles, whether 2 or 
3 parameters are empirical fits, even those based on prop- 
erties (the Dehnen-McLaughlin) of the phase-space density 
profile whose origin is still poorly understood. Therefore the 
only current confidence can be given to those profiles which 
have been fit to the highest resolution simulations and over 
the widest range of halos encountered in N-body simulations. 
Clearly these two criteria are in conflict since simulations 
at the resolution of GHALO are too expensive to allow a 
broader study. Therefore, the results presented here should 
be considered as guides only, whose generality remains to be 
tested. Never the less we can consider economy of parame- 
ters and simplicity of functional form as guiding principles 
in the search for suitable profile functions to describe the 
end state of gravitational collapse. All the profiles we fit 
here (Table 1 and residuals in Figure 2) meet these subjec- 
tive criteria, having at most 3 free parameters and simple 
functional forms. 

While the phase-space density estimated by pa~ 3 is ob- 
served to follow a power law in radius of slope —1.84, its 
meaning is less clear since at small radii it is limited by 
resolution of the estimator and at larger radii it becomes 
dominated by subhalos. Using the more sophisticated En- 
BiD PSD estimator we find that the radial profile is steeper 
with an index of about —2, but that it is not as perfect a 
power-law as seen in pa~ 3 (r). 

As a final comment, we note that in large galaxies, the 
inner structure and shape of the dark matter halo has likely 
been altered over time by the baryons via a range of physical 
effects, including dissipation, energy transfer from sinking 
massive objects, binary black holes, bar-halo interactions, 
turbulent gas motions and more. Simulations that follow 
the baryonic components together with the dark matter will 
resolve these additional questions in the coming years. 
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